This invention relates to stereo imaging velocimetry. Specifically, this invention relates to a system and method for measuring three-dimensional velocities at a plurality of points in a fluid seeded with tracer particles.
Stereo imaging velocimetry seeks to provide a three-dimensional measurement of the velocity of a fluid. A recent example of a stereo imaging velocimetry system is shown in U.S. Pat. No. 5,905,568 of May 18, 1999 which is incorporated by reference herein in its entirety. In this system, a fluid-flow is analyzed using particle tracking. Particle tracking seeks to assign velocity vectors to sequences of still image frames of particles suspended in a transparent fluid, air or gas experiment. As used here, a particle track is a time-sequence of images of the same particle. This requires that high quality particle images be obtained from a system of moving particles. Unfortunately, in practice, images are often contaminated by a variety of noise sources, which introduces significant errors in the velocity measurements. Such errors can cause noise to be misidentified as particles and can cause particles to be misidentified as background noise. A failure to identify a single particle in only one of a sequence of image frames can result in the entire track of the particle being lost.
Examples of noise include stationary noise and moving noise. The primary source of stationary noise may be the container which holds the particles and the experiment. Typically the particles are illuminated with a bright light source, such as a strobe light, which may produce strong glare spots from the container bottom, its bounding side walls, and the liquid surface. Weaker stationary images may also appear within the central field of view due to reflections from overhead lights for example.
Moving noise is generally a reflection of the particles from the surfaces of the container. Moving noise typically occurs at the periphery of the system and as a result may not pose a major problem during processing of image tracks. However stationary noise may introduce image clutter throughout the image which complicates the tracking of tracer particles.
The traditional approach to removing stationary noise under these conditions is to perform time-average, background image subtractions. In this procedure, each image frame is considered part of a shorter fixed-length subsequence of the original time-sequence of image frames, with a xe2x80x98Frame of Interestxe2x80x99 image centered in the subsequence. The xe2x80x98Frame of Interestxe2x80x99 is defined for a given subsequence to be that image frame in the subsequence from which we subtract the time-averaged image calculated using all images in the subsequence. FIG. 1 shows a subsequence 10 which is five image frames in length and whose Frame of Interest 12 is centered in the interval. The image frames comprising the subsequence are averaged together on a pixel-by-pixel basis, and the resultant time-averaged image frame 14 is subtracted from the target image. The net effect to the resulting frame of interest 16 is that rapidly moving particles contribute little to the pixels along their path on average and, therefore are not removed when the average value is subtracted. On the other hand, stationary noise is subtracted out because it is present in every image frame.
Unfortunately when attempting to isolate moving particles by subtracting out a time-average subsequence, two major problems occur. First, requiring that particles be moving relatively rapidly means that slowly-moving particles will be subtracted out as background noise. Second, even rapidly-moving particles will tend to attenuate themselves since they also appear within the subsequence.
The following examples illustrate that these effects can be serious. As shown in FIG. 2, a subsequence of N images 22 beginning at time txcex1, may be taken from a longer sequence of M images 20. The subsequence of intensity values of the pixel 24 located at row xe2x80x98rxe2x80x99 and column xe2x80x98cxe2x80x99 within the subsequence of images 22 may be denoted by Sxcex1(r,c). Since Sxcex1(r,c) is chronologically ordered, the subscript xe2x80x98kxe2x80x99 may be used to denote the kth value of the time sequence: Sxcex1k(r,c). The pixel values may also be placed in an ascending order which may be denoted as Ixcex1(r,c). With this notation, Ixcex1k(r,c) refers to the kth smallest intensity value in Sxcex1(r,c). Also with this notation Ixcex1N(r,c) refers to the largest intensity value.
When a single particle image passes over a specific pixel during the subsequence, the given particle may intersect with this pixel in one or more frames captured during the time period of the subsequence. To simplify the example, it is assumed that no other particle intersects this specific pixel during the time period in question. In addition, gray scale pixel values may range between 0 (black) and 255 (white) and particles are assumed to be white on a black background. However, as shown herein in the drawings, particles are depicted as black spots on a white background to enable the particles in the drawings to be more easily viewed.
In a first example, a subsequence showing a relatively fast moving particle is processed using time averaging. Equation 1 shows an example of the unprocessed ordered subsequence of intensity values (denoted as IF) for a specific pixel of a subsequence.
IF=[13, 45, 67, *126, *244]xe2x80x83xe2x80x83[1]
Here the intensity values marked with an xe2x80x9c*xe2x80x9d correspond to the fast moving particle. The average of the intensity values for IF is 99. Equation 2 shows an example of a processed ordered subsequence of intensity values (denoted as I2F) which corresponds to the unprocessed ordered subsequence IF.
I2F=max{IFxe2x88x92average(IF), 0}=[0, 0, 0, *27, *145]xe2x80x83xe2x80x83[2]
Here the intensity values for I2F are obtained from IF by first subtracting out the time-averaged sequence value, and second setting negative values to zero. In this example the processed subsequence of pixel values I2F, while greatly attenuated, still retain information of the particle passing within each image frame where it was actually present. The fluctuation in the background intensity is a result of low-level noise.
Unfortunately using time-averaged background subtraction for a particle moving much more slowly than the fast moving particle in the previous example may produce undesirable results. In a second example, a subsequence showing a slower moving particle which is not as well lit is processed using time averaging. Such a particle for example may be moving towards the camera along the sides of the experimental chamber. Equation 3 shows an example of the unprocessed ordered subsequence of intensity values (denoted as IS) for a specific pixel of such a subsequence.
IS=[70 ,*120, *143, *168, *176]xe2x80x83xe2x80x83[3]
Here the intensity values marked with an xe2x80x9c*xe2x80x9d correspond to the slow moving particle. The average of the intensity values for IS is 137. Since the particle is moving more slowly, it is present at the pixel under consideration in four of the five frames of the subsequence.
Equation 4 shows an example of a processed ordered subsequence of intensity values (denoted as I2S) which corresponds to the unprocessed ordered subsequence IS.
I2S=max{ISxe2x88x92average(IS), 0}=[0, *0, *8, 33, *41]xe2x80x83xe2x80x83[4]
As in the previous example, the intensity values for I2S are obtained from IS by first subtracting out the time-averaged sequence value, and second setting negative values to zero. In this example the intensity value of *0 in Equation 4 represents a frame when the particle intersected the pixel. However, it has been set to zero (background) due to this prior art method of time-averaged background subtraction.
In addition the use of a threshold for significance would most likely cause the intensify value of xe2x80x9c*8xe2x80x9d in Equation 4 to also be set to zero, further reducing the number of image frames which correctly identify the particle. If the Frame of Interest is one of the frames with a zero or background pixel value, then the particle will disappear from that frame at that pixel location. Unfortunately, when processing frames, the loss of even one frame in the sequence of frames constituting the velocity measurement is enough to cause the loss of that track. When bounded search regions are used to place limits on the size of the search space and thereby reduce the complexity of the tracking algorithm, the opportunity to loose particle tracks may be further increased as a result of time-average background subtraction.
When the velocity measurement is calculated using four frames, there are then four opportunities for a dropped particle to disrupt the identification of the track. The number of tracks in an experiment will also multiply this problem as a result of a phenomenon referred to as track cannibalization. With this phenomenon, not only are the fragmented tracks lost, but also the unassigned particle images which remain behind interact with surrounding tracks during their assignment to further reduce the number of tracks correctly assigned.
In addition time-averaged background subtraction of a subsequence that does not contain a particle may also produce undesirable results. In a third example, a subsequence showing stationary noise of the same magnitude as particle signals is processed using time-averaged background subtraction. Equation 3 shows an example of the unprocessed ordered subsequence of intensity values (denoted as IN) for a specific pixel which experiences only stationery noise during the subsequence of interest.
IN=[0, 13, 45, 51, 67]xe2x80x83xe2x80x83[5]
Here no intensity values are marked with an xe2x80x9c*xe2x80x9d as none of the values corresponds to a moving particle. The average of the intensity values for IN is 35. Subsequence IN includes values which fluctuate quite a bit. This is not uncommon in tracking experiments. For example imperfect synchronization between the strobe and the recording system may be a source of fluctuation of pixel values for otherwise stationary signals. Hence it often common to see background values fluctuate such that intensity values in one or two frames may be well below the intensity values in other frames as shown in Equation 5. In addition glare spots may also be a source of fluctuations in intensity values. Attempts to subtract them out invariably leave xe2x80x9cscintillationsxe2x80x9d throughout the region of the glare.
Equation 6 shows an example of a processed ordered subsequence of intensity values (denoted as I2N) which corresponds to the unprocessed ordered subsequence IN.
I2N=max{INxe2x88x92average(IN), 0}=[0, 0, 10, 16, 32]xe2x80x83xe2x80x83[6]
As in the previous example, the intensity values for I2N are obtained from IN by first subtracting out the time-averaged sequence value, and second setting negative values to zero. In this example the resulting non zero intensity values for the stationary noise are of the same magnitude as the particle signals shown in the previous examples. Any attempt to filter these values out will only server to magnify the problems outlined above concerning loss of tracks. Furthermore, the values depicted in this example are relatively low for stationary noise, which values can run the full range along the boundaries of the containing vessel. The number of lost tracks would be magnified by intensity values for stationary noise which fluctuates over a greater range.
Consequently there exists a need for a system and method of processing image frames in a stereo imaging velocimetry system which is more accurate at identifying true particle image tracks.
In addition there exists a need for a system and method of processing image frames in a stereo imaging velocimetry system which minimizes lost tracks due to slower moving particles.
There further exists a need for a system and method of processing image frames in a stereo imaging velocimetry system which minimizes lost tracks due to stationary noise.
It is an object of an exemplary form of the present invention to provide a system and method for measuring three-dimensional velocities at a plurality of points in a fluid.
It is a further object of an exemplary form of the present invention to provide a system and method for collecting quantitatively, three-dimensional flow data from an optically transparent fluid having tracer particles.
It is a further object of the present invention to provide a stereo imaging velocimetry system and method with improved accuracy.
It is a further object of the present invention to provide a stereo imaging velocimetry system and method which is more accurate at identifying true particle image tracks.
It is a further object of the present invention to provide a stereo imaging velocimetry system and method which minimizes lost tracks due to slower moving particles.
It is a further object of the present invention to provide a stereo imaging velocimetry system and method which minimizes lost tracks due to stationary noise.
Further objects of the present invention will be made apparent in the following Best Modes for Carrying Out Invention and the appended Claims.
The foregoing objects may be accomplished in an exemplary embodiment of the invention by a system that measures three-dimensional velocities at a plurality of points in a fluid seeded with tracer particles. The system may include at least two cameras positioned to view the fluid seeded with tracer particles. The two cameras may be positioned approximately perpendicular to one another to provide two, two-dimensional views of the tracer particles seeded in the fluid. A coordinate system for each camera may be defined based on a world coordinate system. A two-dimensional view of the fluid seeded with the tracer particles may be captured for each camera as a sequence of image frames.
The captured image frames may be processed to remove noise using time-averaged background subtraction with outlier rejection and spike-removal filtering. The system may stereo match particle tracks in the filtered image frames for each camera view and determine three dimensional coordinates of tracer particle tracks. The particle tracking may include using a neural net with simulated annealing to arrive at globally optimal assignments of tracks. The system may further calculate velocities at a plurality of points in the fluid responsive to the three dimensional coordinates of the tracer particle tracks and known time durations between image frames.
An exemplary embodiment of the system may further include calibrating the cameras using a calibration grid. For camera views of the calibration grid that are missing calibration points from the grid, the system may be operative to perform warping transformations which are used to estimate the location of missing points on a calibration grid. The calibration of the cameras may then proceed responsive to both actual and estimated calibration grid points.